Attempts: Portforlio Optimization under semi-G-structure
Yifan Li
04/02/2022
1 Introduction
2 Preparation
source("Port_opt_semi-G-structure-functions.R")##
## Attaching package: 'CVXR'## The following object is masked from 'package:matrixcalc':
##
## vec## The following object is masked from 'package:stats':
##
## power## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'## Also defined by 'Rmpfr'## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.1.1 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::id() masks CVXR::id()
## x purrr::is_vector() masks CVXR::is_vector()
## x dplyr::lag() masks stats::lag()##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutsd.mat0 <- matrix(c(1,3,
2,4), 2, 2, byrow = TRUE)
mu.vec <- matrix(c(1, 2), 2, 1)
rho.vec <- c(0.2,0.4)3 Implementation
3.1 P: min variance
(A1 <- A.f(sd.mat0, rho.vec, par.weight = 1/2, i=1))## [,1] [,2]
## [1,] 5.0 2.6
## [2,] 2.6 10.0(A2 <- A.f(sd.mat0, rho.vec, par.weight = 3/4, i=2))## [,1] [,2]
## [1,] 1.75 0.85
## [2,] 0.85 4.50is.positive.definite(A1)## [1] TRUEis.positive.definite(A2)## [1] TRUE# Variables minimized over
x <- Variable(2)
# Problem definition
objective <- Minimize(quad_form(x, A2/2))
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)
# Problem solution
solution <- solve(prob)
solution$status## [1] "optimal"solution$value## [1] 0.785989solution$getValue(x)## [,1]
## [1,] 0.8021978
## [2,] 0.1978022#average allocation
#achieve smallest variance uncertainty#consider another problem
# Variables minimized over
x <- Variable(2)
t <- 1
# Problem definition
objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)
# Problem solution
solution <- solve(prob)
solution$status## [1] "optimal"solution$value## [1] 0.9102041solution$getValue(x)## [,1]
## [1,] 0.6530612
## [2,] 0.3469388(H.opt <- solution$getValue(quad_form(x, A1)))## [1] 4.514286(mu.opt <- solution$getValue(t(mu.vec)%*%x))## [1] 1.3469393.2 Q: efficient frontier
3.2.1 A1
#draw the efficient frontier
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
# x.opt.mat <- matrix(NA,
# nrow = length(t.seq), ncol=2)
for (i in seq_along(t.seq)){
t <- t.seq[i]
x <- Variable(2)
# Problem definition
objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)
# Problem solution
solution <- solve(prob)
#x.opt.mat[i,] <- as.numeric(solution$getValue(x))
H.opt.seq[i] <- solution$getValue(quad_form(x, A1))
mu.opt.seq[i] <- solution$getValue(t(mu.vec)%*%x)
}#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l",
main = "A1 opt problem")
3.2.2 A2
#draw the efficient frontier
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
for (i in seq_along(t.seq)){
t <- t.seq[i]
x <- Variable(2)
# Problem definition
objective <- Minimize(quad_form(x, A2/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)
# Problem solution
solution <- solve(prob)
H.opt.seq[i] <- solution$getValue(quad_form(x, A1))
mu.opt.seq[i] <- solution$getValue(t(mu.vec)%*%x)
}#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l",
main = "A2 opt problem")
EF does depend on the measure we choose and also the tuning parameter (which depends on the preference of the user on the variance uncertainty).
3.3 EF: change par
3.3.1 A1
t.seq <- seq(-10, 10, .2)
a.seq <- seq(0.1, 0.9,.1)
H.opt.mat <- mu.opt.mat <- matrix(NA,
ncol = length(a.seq),
nrow = length(t.seq))
#we can run parallel computing for this part
for (j in seq_along(a.seq)){
a <- a.seq[j]
A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
t <- t.seq[i]
x <- Variable(2)
# Problem definition
objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)
# Problem solution
solution <- solve(prob)
H.opt.mat[i, j] <- solution$getValue(quad_form(x, A1))
mu.opt.mat[i, j] <- solution$getValue(t(mu.vec)%*%x)
}
}
#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat), max(H.opt.mat))
mu.lim <- c(min(mu.opt.mat), max(mu.opt.mat))
for (j in seq_along(a.seq)){
if (j==1){
plot(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
type = "l", xlab = "H.opt", ylab = "mu.opt",
xlim = H.lim, ylim = mu.lim,
main = "Changing tuning par for P2.1")
} else {
lines(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
type = "l")
}
}3.3.2 A2
t.seq <- seq(-10, 10, .2)
a.seq <- seq(2/3, 0.99, length.out = 9)
H.opt.mat2 <- mu.opt.mat2 <- matrix(NA,
ncol = length(a.seq),
nrow = length(t.seq))
#we can run parallel computing for this part
#seq_along(a.seq)
for (j in seq_along(a.seq)){
a <- a.seq[j]
A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
t <- t.seq[i]
x <- Variable(2)
# Problem definition
objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)
# Problem solution
solution <- solve(prob)
H.opt.mat2[i, j] <- solution$getValue(quad_form(x, A1))
mu.opt.mat2[i, j] <- solution$getValue(t(mu.vec)%*%x)
}
}
#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat2), max(H.opt.mat2))
mu.lim <- c(min(mu.opt.mat2), max(mu.opt.mat2))
for (j in seq_along(a.seq)){
if (j==1){
plot(H.opt.mat2[,j],mu.opt.mat2[,j], col=j,lty=j,
type = "l", xlab = "H.opt", ylab = "mu.opt",
xlim = H.lim, ylim = mu.lim,
main = "Changing tuning par for P2.2")
} else {
lines(H.opt.mat2[,j],mu.opt.mat2[,j], col=j,lty=j,
type = "l")
}
}4 Extended cases
4.1 Objective function
#consider the rho as negative
#consider this problem from theoretical level#check whether the optimal weight depend on the tuning parameter#plot the range of the functions
#H is the objective function
#draw a 3D surface plot
sig1 <- seq(0.1, 4, .01)
sig2 <- seq(2, 5, .01)
z <- outer(sig1, sig2, H)
contour(sig1, sig2, z, xlab="sig1", ylab="sig2", xlim = range(sig1), ylim = range(sig2))
#abline(a=0, b=0.5, col=2, xlim = range(sig1))
#abline(a=0, b=1/0.5, col=2, ylim = range(sig2))which.min(z)## [1] 16which.max(z)## [1] 391plot(function(x) H(x, sig2=2), from = min(sig1), to = max(sig1),
xlab = "sig1")
plot(function(x) H(sig1=2, sig2=x), from = min(sig2), to = max(sig2),
xlab = "sig2")
x <- sig1
y <- sig2
fig <- plot_ly(
type = 'surface',
contours = list(
z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
x = ~x,
y = ~y,
z = ~z)
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 20),
zaxis = list(nticks = 20),
camera = list(eye = list(x = 0, y = -1, z = 0.5)),
aspectratio = list(x = .9, y = .8, z = .6)))
fig#our current guess is, for
#check different cases
#the objective function may no longer be a convex one#use range.set
range.set(x1=0.2, print.ind = TRUE, rho.intvl = c(-0.8, -0.4))## i = 1 , sig1 = 0.5 , sig2 = 1 , rho = -0.8 , H = 0.522
## i = 2 , sig1 = 0.5 , sig2 = 2 , rho = -0.8 , H = 2.314
## i = 3 , sig1 = 1 , sig2 = 1 , rho = -0.8 , H = 0.424
## i = 4 , sig1 = 1 , sig2 = 2 , rho = -0.8 , H = 2.088
## i = 5 , sig1 = 0.5 , sig2 = 1 , rho = -0.4 , H = 0.586
## i = 6 , sig1 = 0.5 , sig2 = 2 , rho = -0.4 , H = 2.442
## i = 7 , sig1 = 1 , sig2 = 1 , rho = -0.4 , H = 0.552
## i = 8 , sig1 = 1 , sig2 = 2 , rho = -0.4 , H = 2.344## [1] 0.424 2.442range.set(x1=0.5, print.ind = TRUE, rho.intvl = c(-0.8, -0.4))## i = 1 , sig1 = 0.5 , sig2 = 1 , rho = -0.8 , H = 0.1125
## i = 2 , sig1 = 0.5 , sig2 = 2 , rho = -0.8 , H = 0.6625
## i = 3 , sig1 = 1 , sig2 = 1 , rho = -0.8 , H = 0.1
## i = 4 , sig1 = 1 , sig2 = 2 , rho = -0.8 , H = 0.45
## i = 5 , sig1 = 0.8 , sig2 = 1 , rho = -0.8 , H = 0.09
## i = 6 , sig1 = 0.8 , sig2 = 2 , rho = -0.8 , H = 0.52
## i = 7 , sig1 = 0.5 , sig2 = 1 , rho = -0.4 , H = 0.2125
## i = 8 , sig1 = 0.5 , sig2 = 2 , rho = -0.4 , H = 0.8625
## i = 9 , sig1 = 1 , sig2 = 1 , rho = -0.4 , H = 0.3
## i = 10 , sig1 = 1 , sig2 = 2 , rho = -0.4 , H = 0.85
## i = 11 , sig1 = 0.8 , sig2 = 1 , rho = -0.4 , H = 0.25
## i = 12 , sig1 = 0.8 , sig2 = 2 , rho = -0.4 , H = 0.84
## i = 13 , sig1 = 0.8 , sig2 = 1 , rho = -0.4 , H = 0.25
## i = 14 , sig1 = 0.8 , sig2 = 2 , rho = -0.4 , H = 0.84## [1] 0.0900 0.8625range.set(x1=0.8, print.ind = TRUE, rho.intvl = c(-0.8, -0.4))## i = 1 , sig1 = 0.5 , sig2 = 1 , rho = -0.8 , H = 0.072
## i = 2 , sig1 = 0.5 , sig2 = 2 , rho = -0.8 , H = 0.064
## i = 3 , sig1 = 0.5 , sig2 = 1.6 , rho = -0.8 , H = 0.0576
## i = 4 , sig1 = 1 , sig2 = 1 , rho = -0.8 , H = 0.424
## i = 5 , sig1 = 1 , sig2 = 2 , rho = -0.8 , H = 0.288
## i = 6 , sig1 = 1 , sig2 = 1.6 , rho = -0.8 , H = 0.3328
## i = 7 , sig1 = 0.5 , sig2 = 1 , rho = -0.4 , H = 0.136
## i = 8 , sig1 = 0.5 , sig2 = 2 , rho = -0.4 , H = 0.192
## i = 9 , sig1 = 0.5 , sig2 = 1.6 , rho = -0.4 , H = 0.16
## i = 10 , sig1 = 0.5 , sig2 = 1.6 , rho = -0.4 , H = 0.16
## i = 11 , sig1 = 1 , sig2 = 1 , rho = -0.4 , H = 0.552
## i = 12 , sig1 = 1 , sig2 = 2 , rho = -0.4 , H = 0.544
## i = 13 , sig1 = 1 , sig2 = 1.6 , rho = -0.4 , H = 0.5376
## i = 14 , sig1 = 1 , sig2 = 1.6 , rho = -0.4 , H = 0.5376## [1] 0.0576 0.5520#plot the maximum variance
x1.seq <- seq(-1, 2, .01)
#x1.seq <- seq(0, 1, .01)
rho.intvl1 <- c(-0.8,-0.4)
h.seq <- sapply(x1.seq, function(x){
re <- range.set(x1=x, rho.intvl = rho.intvl1)
re[2]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var",
main = paste0("rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
#the max variance
#it is still a convex function
#not only for positive x1rho.intvl1 <- c(-0.8,0.4)
h.seq <- sapply(x1.seq, function(x){
re <- range.set(x1=x, rho.intvl = rho.intvl1)
re[2]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var",
main = paste0("rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
rho.intvl1 <- c(-0.8,0.4)
h.seq <- sapply(x1.seq, function(x){
re <- range.set(x1=x, rho.intvl = rho.intvl1)
re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "min var",
main = paste0("rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
rho.intvl1 <- c(-0.2,0.2)
h.seq <- sapply(x1.seq, function(x){
re <- range.set(x1=x, rho.intvl = rho.intvl1)
re[2]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var",
main = paste0("rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
#plot the min variance
x1.seq <- seq(-1, 2, .01)
#x1.seq <- seq(0, 1, .01)
rho.intvl1 <- c(-0.8,-0.4)
h.seq <- sapply(x1.seq, function(x){
re <- range.set(x1=x, rho.intvl = rho.intvl1)
re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "min var",
main = paste0("rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
#the min variance
#it is not a convex function
#not only for positive x1x1.seq <- seq(0.9, 1.1, .01)
rho.intvl1 <- c(-0.8,-0.4)
h.seq <- sapply(x1.seq, function(x){
re <- range.set(x1=x, rho.intvl = rho.intvl1)
re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "min var",
main = paste0("rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
x1.seq <- seq(-1, 2, .01)
rho.intvl1 <- c(-0.8, -0.5)
h.seq <- sapply(x1.seq, function(x){
re <- range.set(x1=x, rho.intvl = rho.intvl1)
re[2]-re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var",
main = paste0("rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
#It may not be convex
#but it should still have a global minimum4.2 P: min variance
4.2.1 A1
#run the numerical optimization for this function
#so far we do not have a closed form solution
#direct minimization of H w.r.t x1
#both report the variance and current mean
E1 <- function(x, lam = 0.5, rho.intvl = c(-0.8, -0.5)){
re <- range.set(x1=x, rho.intvl = rho.intvl)
lam*re[2]+(1-lam)*re[1]
}
(opt.re <- optimize(E1, interval = c(-50,50)))## $minimum
## [1] 0.7071823
##
## $objective
## [1] 0.2370166x1.opt <- opt.re$minimum
(x.opt <- c(x1.opt, 1-x1.opt))## [1] 0.7071823 0.2928177(obj.opt <- opt.re$objective)## [1] 0.23701664.2.2 A2
E2 <- function(x, kap = 0.5, rho.intvl = c(-0.8, -0.5)){
re <- range.set(x1=x, rho.intvl = rho.intvl)
cen <- (re[2]+re[1])/2
amb <- re[2]-re[1]
kap*cen+(1-kap)*amb
}
(opt.re <- optimize(E2, c(-50, 50)))## $minimum
## [1] 0.7173543
##
## $objective
## [1] 0.3098996x1.opt <- opt.re$minimum
(x.opt <- c(x1.opt, 1-x1.opt))## [1] 0.7173543 0.2826457(obj.opt <- opt.re$objective)## [1] 0.30989964.3 Q: efficient frontier
- change the objective function
- subtract by \(tx\mu\)
- then also draw the EF
4.3.1 A1
#draw the efficient frontier
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
mu.vec <- as.numeric(mu.vec)
# x.opt.mat <- matrix(NA,
# nrow = length(t.seq), ncol=2)
rho.intvl1 <- c(-0.8, -0.5)
#objective function
E1 <- function(x, lam = 0.5, rho.intvl = rho.intvl1){
re <- range.set(x1=x, rho.intvl = rho.intvl)
lam*re[2]+(1-lam)*re[1]
}
for (i in seq_along(t.seq)){
t <- t.seq[i]
obj <- function(x) E1(x) - t*sum(mu.vec*x)
opt.re <- optimize(obj, interval = c(-50,50))
x1.opt <- opt.re$minimum
x.opt <- c(x1.opt, 1-x1.opt)
H.opt.seq[i] <- E1(x1.opt)
mu.opt.seq[i] <- sum(mu.vec*x.opt)
}
#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l",
main = paste0("Q1 opt problem with ", "rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
4.3.2 A2
#draw the efficient frontier
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
mu.vec <- as.numeric(mu.vec)
# x.opt.mat <- matrix(NA,
# nrow = length(t.seq), ncol=2)
#objective function
rho.intvl1 <- c(-0.8, -0.5)
E2 <- function(x, kap = 0.5, rho.intvl = rho.intvl1){
re <- range.set(x1=x, rho.intvl = rho.intvl)
cen <- (re[2]+re[1])/2
amb <- re[2]-re[1]
kap*cen+(1-kap)*amb
}
for (i in seq_along(t.seq)){
t <- t.seq[i]
obj <- function(x) E2(x) - t*sum(mu.vec*x)
opt.re <- optimize(obj, interval = c(-50,50))
x1.opt <- opt.re$minimum
x.opt <- c(x1.opt, 1-x1.opt)
H.opt.seq[i] <- E2(x1.opt) #to change
mu.opt.seq[i] <- sum(mu.vec*x.opt)
}
#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l",
main = paste0("Q1 opt problem with ", "rho.intvl = [", rho.intvl1[1],
",", rho.intvl1[2], "]"))
4.4 EF: change par
t.seq <- seq(-5, 5, .1)
a.seq <- seq(0.1, 0.9,.1)4.4.1 A1
H.opt.mat <- mu.opt.mat <- matrix(NA,
ncol = length(a.seq),
nrow = length(t.seq))
#we can run parallel computing for this part
for (j in seq_along(a.seq)){
a <- a.seq[j]
E1 <- function(x, lam = a, rho.intvl = c(-0.8, 0.5)){
re <- range.set(x1=x, rho.intvl = rho.intvl)
lam*re[2]+(1-lam)*re[1]
}
#A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
t <- t.seq[i]
obj <- function(x) E1(x) - t*sum(mu.vec*x)
opt.re <- optimize(obj, interval = c(-50,50))
x1.opt <- opt.re$minimum
x.opt <- c(x1.opt, 1-x1.opt)
H.opt.mat[i, j] <- E1(x1.opt)
mu.opt.mat[i, j] <- sum(mu.vec*x.opt)
}
}## j = 1
## j = 2
## j = 3
## j = 4
## j = 5
## j = 6
## j = 7
## j = 8
## j = 9#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat), max(H.opt.mat))
mu.lim <- c(min(mu.opt.mat), max(mu.opt.mat))
for (j in seq_along(a.seq)){
if (j==1){
plot(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
type = "l", xlab = "H.opt", ylab = "mu.opt",
xlim = H.lim, ylim = mu.lim,
main = "Changing tuning par for Q1")
} else {
lines(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
type = "l")
}
}
4.4.2 A2
H.opt.mat <- mu.opt.mat <- matrix(NA,
ncol = length(a.seq),
nrow = length(t.seq))
#we can run parallel computing for this part
for (j in seq_along(a.seq)){
a <- a.seq[j]
E2 <- function(x, kap = a, rho.intvl = c(-0.8, 0.5)){
re <- range.set(x1=x, rho.intvl = rho.intvl)
cen <- (re[2]+re[1])/2
amb <- re[2]-re[1]
kap*cen+(1-kap)*amb
}
#A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
t <- t.seq[i]
obj <- function(x) E2(x) - t*sum(mu.vec*x)
opt.re <- optimize(obj, interval = c(-50,50))
x1.opt <- opt.re$minimum
x.opt <- c(x1.opt, 1-x1.opt)
H.opt.mat[i, j] <- E1(x1.opt)
mu.opt.mat[i, j] <- sum(mu.vec*x.opt)
}
}## j = 1
## j = 2
## j = 3
## j = 4
## j = 5
## j = 6
## j = 7
## j = 8
## j = 9#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat), max(H.opt.mat))
mu.lim <- c(min(mu.opt.mat), max(mu.opt.mat))
for (j in seq_along(a.seq)){
if (j==1){
plot(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
type = "l", xlab = "H.opt", ylab = "mu.opt",
xlim = H.lim, ylim = mu.lim,
main = "Changing tuning par for Q1")
} else {
lines(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
type = "l")
}
}
4.5 How does the optimal port depend on the paramter?
4.5.1 A1
4.5.1.1 sig1.intvl
#x1 vs sig1.intvl
x.seq <- sig1.low.seq <- seq(0.5, 4, .1)
y.seq <- sig1.amb.seq <- seq(0.1, 2, .1)
S <- function(a,b){
#both report the variance and current mean
E1 <- function(x, lam = 0.5, rho.intvl = c(0.5, 0.5),
sig1.intvl = lowamb2intvl(a,b),
sig2.intvl = c(2, 2)){
re <- range.set(x1=x, rho.intvl = rho.intvl,
sig1.intvl = sig1.intvl,
sig2.intvl = sig2.intvl)
lam*re[2]+(1-lam)*re[1]
}
opt.re <- optimize(E1, interval = c(-50,50))
x1.opt <- opt.re$minimum
return(x1.opt)
}
z <- matrix(NA, nrow = length(x.seq), ncol = length(y.seq))
for (i in seq_along(x.seq)){
for (j in seq_along(y.seq)){
z[i,j] <- S(x.seq[i], y.seq[j])
}
}
#z <- outer(x.seq, y.seq, S)
#fig1 <- plot3d(x.seq, y.seq, z.mat)
#rglwidget(elementId = "plot3drgl")
x <- x.seq
y <- y.seq
fig <- plot_ly(
type = 'surface',
# contours = list(
# z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
x = ~x,
y = ~y,
z = ~z)
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 10, title = "sig1.low"),
yaxis = list(nticks = 10, title = "sig1.amb"),
zaxis = list(nticks = 10, title = "x1.opt"),
camera = list(eye = list(x = -1, y = -1, z = 0)),
aspectratio = list(x = .9, y = .8, z = 0.6)))
figAs \(\underline{\sigma}_1\) increases or the ambiguity of variance (of asset one) increases, the proportion on \(x1\) will decrease.
#x1 vs sig1.intvl
x.seq <- sig1.low.seq <- seq(0.5, 4, .1)
y.seq <- sig1.amb.seq <- seq(0.1, 2, .1)
S <- function(a,b){
#both report the variance and current mean
E1 <- function(x, lam = 0.5, rho.intvl = c(-0.5, -0.5),
sig1.intvl = lowamb2intvl(a,b),
sig2.intvl = c(2, 2)){
re <- range.set(x1=x, rho.intvl = rho.intvl,
sig1.intvl = sig1.intvl,
sig2.intvl = sig2.intvl)
lam*re[2]+(1-lam)*re[1]
}
opt.re <- optimize(E1, interval = c(-50,50))
x1.opt <- opt.re$minimum
return(x1.opt)
}
z <- matrix(NA, nrow = length(x.seq), ncol = length(y.seq))
for (i in seq_along(x.seq)){
for (j in seq_along(y.seq)){
z[i,j] <- S(x.seq[i], y.seq[j])
}
}
#z <- outer(x.seq, y.seq, S)
#fig1 <- plot3d(x.seq, y.seq, z.mat)
#rglwidget(elementId = "plot3drgl")
x <- x.seq
y <- y.seq
fig <- plot_ly(
type = 'surface',
# contours = list(
# z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
x = ~x,
y = ~y,
z = ~z)
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 10, title = "sig1.low"),
yaxis = list(nticks = 10, title = "sig1.amb"),
zaxis = list(nticks = 10, title = "x1.opt"),
camera = list(eye = list(x = -1, y = 0, z = 0)),
aspectratio = list(x = .9, y = .8, z = .6)))
fig#x1 vs sig1.intvl
x.seq <- sig1.low.seq <- seq(0.5, 4, .1)
y.seq <- sig1.amb.seq <- seq(0.1, 2, .1)
S <- function(a,b){
#both report the variance and current mean
E1 <- function(x, lam = 0.5, rho.intvl = c(-0.5, 0.5),
sig1.intvl = lowamb2intvl(a,b),
sig2.intvl = c(2, 2)){
re <- range.set(x1=x, rho.intvl = rho.intvl,
sig1.intvl = sig1.intvl,
sig2.intvl = sig2.intvl)
lam*re[2]+(1-lam)*re[1]
}
opt.re <- optimize(E1, interval = c(-50,50))
x1.opt <- opt.re$minimum
return(x1.opt)
}
z <- matrix(NA, nrow = length(x.seq), ncol = length(y.seq))
for (i in seq_along(x.seq)){
for (j in seq_along(y.seq)){
z[i,j] <- S(x.seq[i], y.seq[j])
}
}
#z <- outer(x.seq, y.seq, S)
#fig1 <- plot3d(x.seq, y.seq, z.mat)
#rglwidget(elementId = "plot3drgl")
x <- x.seq
y <- y.seq
fig <- plot_ly(
type = 'surface',
# contours = list(
# z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
x = ~x,
y = ~y,
z = ~z)
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 10, title = "sig1.low"),
yaxis = list(nticks = 10, title = "sig1.amb"),
zaxis = list(nticks = 10, title = "x1.opt"),
camera = list(eye = list(x = -1, y = 0, z = 0)),
aspectratio = list(x = .9, y = .8, z = .6)))
fig4.5.1.2 rho.intvl
#x1 vs sig1.intvl
x.seq <- rho.l.seq <- seq(-0.9, 0.9, .01)
y.seq <- rho.r.seq <- seq(-0.9, 0.9, .01)
load("z_rho.Rdata")S <- function(a,b){
if (a>b) {
return(NA)
} else {
E1 <- function(x, lam = 0.5, rho.intvl = c(a,b),
sig1.intvl = c(0.5,1),
sig2.intvl = c(1,2)){
re <- range.set(x1=x, rho.intvl = rho.intvl,
sig1.intvl = sig1.intvl,
sig2.intvl = sig2.intvl)
lam*re[2]+(1-lam)*re[1]
}
opt.re <- optimize(E1, interval = c(-50,50))
x1.opt <- opt.re$minimum
return(x1.opt)
}
}
z <- matrix(NA, ncol = length(x.seq), nrow = length(y.seq))
for (j in seq_along(x.seq)){
for (i in seq_along(y.seq)){
z[i,j] <- S(x.seq[j], y.seq[i])
}
}
#z <- outer(x.seq, y.seq, S)
#save(z, file = "z_rho.Rdata")x <- x.seq
y <- y.seq
fig <- plot_ly(
type = 'surface',
# contours = list(
# z = list(show = TRUE, start = 0, end = max(z), size = 0.02)),
x = ~x,
y = ~y,
z = ~z
)
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 10, title = "rho.low"),
yaxis = list(nticks = 10, title = "rho.up"),
zaxis = list(nticks = 10, title = "x1.opt"),
camera = list(eye = list(x = -1, y = 0, z = 0)),
aspectratio = list(x = .9, y = .8, z = .6)))
fig4.5.1.3 tuning parameter
(learn how to better interpret the results.)
5 Eigenvalue decomposition
(show the geometric result)
6 Directly visualize the function
6.1 Setup
#check sigma
d <- 1e-4
x.seq <- seq(d, 1-d,d)
# sd1 <- c(0.25, 1)
# sd2 <- c(0.75, 2)
sd1 <- c(0.5, 2)
sd2 <- c(2.5, 3)
sd1.c <- mean(sd1)
sd2.c <- mean(sd2)
s.set <- rep(sd1, each=2)/rep(sd2, 2)
rhol <- -0.8
sd1[1]*sd2[2] < sd1[2]*sd2[1]## [1] TRUErhol^2 > (sd1[1]*sd2[1])/(sd1[2]*sd2[2])## [1] TRUErhol^2 < sd1[1]/sd1[2]## [1] FALSE# s.set.sort <- sort(c(s.set*(-rhol),
# 1/s.set*(-rhol)
# #s.set/(-rhol),
# #1/(s.set*(-rhol))
# ))
# r.thres <- s.set.sort
# x.thres <- 1/(1+r.thres)
# r1 <- sd2[1]/(sd1[2]*(-rhol))
# r2 <- (sd1[1]*(-rhol))/sd2[2]#important note
#the meaning of lower, upper in optim() is important!
re.mat <- matrix(NA, nrow = length(x.seq), ncol = 2)
val.seq <- numeric(length(x.seq))
for (i in seq_along(x.seq)){
x1 <- x.seq[i]
re <- optim(c(sd1.c,sd2.c), function(s) H(s[1],s[2], rho = rhol, x1 = x1),
lower = c(sd1[1], sd2[1]),
upper = c(sd1[2], sd2[2]),
method = "L-BFGS-B")
#not lower=sd1, upper=sd2
re.mat[i,]<- re$par
#get the min value
val.seq[i] <- re$value
}
#any(is.nan(re.mat[,2]))
#re.mat
summary(re.mat[,1])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5000 0.6668 2.0000 1.4400 2.0000 2.0000summary(re.mat[,2])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.500 2.500 2.500 2.564 2.500 3.000# optim(c(1,1,1), function(x) sum(x^2), lower = c(0,0,0), upper = c(2,2,2), method = "L-BFGS-B")#r1 <- 1/abs(rhol)
#x.thres1 <- c(x.thres, 1/(1+r1))
matplot(x.seq, re.mat, type = "l", main = paste0("min cov ","\n", "with sd1 = ", "[", sd1[1], ",", sd1[2],"], ",
"sd2 = ", "[", sd2[1], ",", sd2[2],"]"), lty = c(1,1))
abline(v = x.thres.f(s.set*(-rhol))[c(1,2)], lty = 2, col = "purple")
abline(v = x.thres.f(s.set/(-rhol))[c(1,3)], lty = 2, col = "blue")
# abline(v = x.thres.f((-rhol)/s.set), lty = 2, col = "brown")
# abline(v = x.thres.f(1/(s.set*(-rhol))), lty = 2, col = "pink")
#abline(v = x.thres.f(s.set/(-rhol)), lty = 1:4, col = 1:4)
#I have got it!plot(x.seq, val.seq, type = "l")
abline(v = x.thres.f(s.set*(-rhol))[c(1,2)], lty = 2, col = "purple")
abline(v = x.thres.f(s.set/(-rhol))[c(1,3)], lty = 2, col = "blue")
t.break <- c(x.thres.f(s.set/(-rhol))[c(3,1)],
x.thres.f(s.set*(-rhol))[c(1,2)])
re.mat.c <- matrix(NA, ncol=2, nrow = length(x.seq))
sdl1 <- sd1[1]
sdr1 <- sd1[2]
sdl2 <- sd2[1]
sdr2 <- sd2[2]
q <- 1/(-rhol)
#check the function
for (i in seq_along(x.seq)){
x1 <- x.seq[i]
x2 <- 1-x1
r <- x2/x1
if (x1 <= t.break[1]){
re.mat.c[i,] <- c(sdr1, sdl2)
} else if (x1 <= t.break[2]) {
re.mat.c[i,] <- c(r*sdl2/q, sdl2)
} else if (x1 <= t.break[3]) {
re.mat.c[i,] <- c(sdl1, sdl2)
} else if (x1 <= t.break[4]) {
re.mat.c[i,] <- c(sdl1, sdl1/(r*q))
} else {
re.mat.c[i,] <- c(sdl1, sdr2)
}
#as x1 increases, r decreases
# if (r <= r.break[1]){
# re.mat.c[i,] <- c(sdl1, sdr2)
# } else if (r <= r.break[2]) {
# re.mat.c[i,] <- c(sdl1, q*sdl1/r)
# } else if (r <= r.break[3]) {
# re.mat.c[i,] <- c(sdl1, sdl2)
# } else if (r <= r.break[4]) {
# re.mat.c[i,] <- c(r*sdl2/q, sdl2)
# } else {
# re.mat.c[i,] <- c(sdr1, sdl2)
# }
}
matplot(x.seq, re.mat.c, type = "l", main = paste0("min cov.check ","\n", "with sd1 = ", "[", sd1[1], ",", sd1[2],"], ",
"sd2 = ", "[", sd2[1], ",", sd2[2],"]"), lty = c(1,1), ylim = c(0.5, 3))
abline(v = x.thres.f(s.set*(-rhol))[c(1,2)], lty = 2, col = "purple")
abline(v = x.thres.f(s.set/(-rhol))[c(1,3)], lty = 2, col = "blue")
r.seq <- (1-x.seq)/x.seq
matplot(log(r.seq), re.mat, type = "l", lty = c(1,1))
re.mat.d <- apply(re.mat, 2, diff)
matplot(x.seq[-1], re.mat.d, type = "l")
ind1 <- which(re.mat.d[,1]<0)
x.t1 <- c(min(x.seq[ind1]), max(x.seq[ind1]))
ind2 <- which(re.mat.d[,2]>0)
x.t2 <- c(min(x.seq[ind2]), max(x.seq[ind2]))x.t <- c(x.t1, x.t2)
x.t## [1] 0.5000 0.7999 0.8620 0.8823#x.t*(-rhol)
r.t <- (1-x.t)/x.t
r.t## [1] 1.0000000 0.2501563 0.1600928 0.1334013r.t*(-rhol)## [1] 0.8000000 0.2001250 0.1280742 0.1067211r.t/(-rhol)## [1] 1.2500000 0.3126953 0.2001160 0.1667517# x.seq <- seq(.01, .99,.001)
re.mat <- matrix(NA, nrow = length(x.seq), ncol = 2)
for (i in seq_along(x.seq)){
x1 <- x.seq[i]
re <- optim(c(1,2), function(s) -H(s[1],s[2], rho = rhol, x1 = x1),
lower = c(sd1[1], sd2[1]),
upper = c(sd1[2], sd2[2]),
method = "L-BFGS-B")
re.mat[i,]<- re$par
}
#re.mat
matplot(x.seq, re.mat, type = "l", main = "max cov", lty = c(1,1))
#r.seq <- (1-x.seq)/x.seq
#matplot(r.seq, re.mat, type = "l")
#re.mat.d <- apply(re.mat, 2, diff)
#matplot(x.seq[-1], re.mat.d, type = "l")tan.theta <- function(x1, rho = -0.5){
x2 <- 1-x1
a <- 2*x1*x2*rho
D <- sqrt((x1^2-x2^2)^2 + a^2)
b <- x1^2-x2^2 + D
a/b
}
plot(tan.theta, from = 0, to = 1)
#tan.theta(0)
tan.theta(1)## [1] 0theta.f <- function(x1, rho = -0.5){
if(x1==0) return(pi/2)
if(x1==1) return(0)
x2 <- 1-x1
a <- 2*x1*x2*rho
D <- sqrt((x1^2-x2^2)^2 + a^2)
b <- x1^2-x2^2 + D
atan(a/b)+pi
}
x <- seq(0,1,.01)
#sapply(x, theta.f)
plot(x, sapply(x, theta.f), type = "l")
#plot(theta.f, from = 0, to = 1)
theta.f(0.99999)## [1] 3.141588# px1 = 0.999
# x2 <- 1-x1
# a <- 2*x1*x2*rho
# D <- sqrt((x1^2-x2^2)^2 + a^2)
# b <- x1^2-x2^2 + D
# atan(a/b) 6.2 Cen-amb graph
Center-ambiguity graph (to see the center and ambiguity of variance interval for a portfolio under different weights)
#check sigma
d <- 1e-4
x.seq <- seq(d, 1-d,d)
sd1 <- c(0.25, 1.5)
sd2 <- c(1, 2)
# sd1 <- c(0.5, 2)
# sd2 <- c(2.5, 3)
sd1.c <- mean(sd1)
sd2.c <- mean(sd2)
s.set <- rep(sd1, each=2)/rep(sd2, 2)
rhol <- -0.8
rhor <- -0.2
sd1[1]*sd2[2] < sd1[2]*sd2[1]## [1] TRUErhol^2 > (sd1[1]*sd2[1])/(sd1[2]*sd2[2])## [1] TRUErhol^2 < sd1[1]/sd1[2]## [1] FALSEtitle.main <- paste0("sd1 = ", "[", sd1[1], ",", sd1[2],"], ", "sd2 = ", "[", sd2[1], ",", sd2[2],"]", ", ",
"rho = ", "[", rhol, ",", rhor,"]")
print(title.main)## [1] "sd1 = [0.25,1.5], sd2 = [1,2], rho = [-0.8,-0.2]"#important note
#the meaning of lower, upper in optim() is important!
# re.mat <- matrix(NA, nrow = length(x.seq), ncol = 2)
val.min.seq <- val.max.seq <- numeric(length(x.seq))
for (i in seq_along(x.seq)){
x1 <- x.seq[i]
#min var
re.min <- optim(c(sd1.c,sd2.c), function(s) H(s[1],s[2], rho = rhol, x1 = x1),
lower = c(sd1[1], sd2[1]),
upper = c(sd1[2], sd2[2]),
method = "L-BFGS-B")
val.min.seq[i] <- re.min$value
#max var
re.max <- optim(c(sd1.c,sd2.c), function(s) -H(s[1],s[2], rho = rhor, x1 = x1),
lower = c(sd1[1], sd2[1]),
upper = c(sd1[2], sd2[2]),
method = "L-BFGS-B")
val.max.seq[i] <- -re.max$value
}
matplot(x.seq, cbind(val.min.seq,val.max.seq), type = "l", lty = c(1,1), main = title.main)
#center
var.cen.seq <- (val.min.seq+val.max.seq)/2
plot(x.seq, var.cen.seq, type = "l",
main = title.main)
#amb
var.amb.seq <- val.max.seq - val.min.seq
plot(x.seq, var.amb.seq, type = "l",
main = title.main)
#plot it
plot(var.cen.seq, var.amb.seq,
main = title.main, type = "p")
#minimizer is very close to each other
matplot(x.seq, cbind(var.cen.seq,var.amb.seq), type = "l", lty = c(1,1), main = title.main)
abline(v = x.seq[which.min(var.cen.seq)], col = 1, lty=2)
abline(v = x.seq[which.min(var.amb.seq)], col = 2, lty=2)
library(ggplot2)
x1 <- x.seq
dat1 <- data.frame(x=var.cen.seq, y=var.amb.seq)
ggplot(dat1, aes(x=x, y=y)) + geom_point(aes(colour = x1)) + xlab("Var.center") + ylab("Var.amb") + ggtitle(title.main) +
scale_colour_gradient(low = "#00FFFF", high = "#000099")
#scale_colour_gradient(low = "lightblue", high = "darkblue")
#use different color to show the value of x1
#similar to efficient frontier It seems that there is a unique minimizer for smallest cen and ambiguity for a large class of objective function (linear combinations of center and ambiguity, use sum of square of two components). The minimizer does not depend much on the format of the function.
The mean part will also change monotonically w.r.t \(x1\), because \[ \mu = x_1\mu_1 + x_2\mu_2 \]
#write a function
#plot.ambcen
plot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, -0.2))
plot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, 0.2))
plot.ambcen.var(sd1 = c(0.5, 2.26), sd2 = c(1,2.03), rho = c(-0.8, -0.4))
#we can notice some "discontinuity"
#it should only come from some numerical issues, and the current x1 may go through some threshold. # plot.ambcen.var(sd1 = c(0.2, 1), sd2 = c(1.2,2.4),
# rho = c(-0.8, -0.2))(Do we also want to consider the mean part in this picture?)
#create a shiny app if needed
#We can also change how the amb-cen curve changes depending on the parameter setup
#to let user change the parameter in a convenient way
#show the optimal x1 by drawing a straight line or using a circle (to better show the objective function)
library(shiny)
ui <- basicPage(
plotOutput("plot1", click = "plot_click"),
verbatimTextOutput("info")
)
server <- function(input, output) {
output$plot1 <- renderPlot({
plot(mtcars$wt, mtcars$mpg)
})
output$info <- renderText({
paste0("x=", input$plot_click$x, "\ny=", input$plot_click$y)
})
}
shinyApp(ui, server)6.3 Min-max var plot
x1 <- x.seq
dat1 <- data.frame(x=val.min.seq, y=val.max.seq)
ggplot(dat1, aes(x=x, y=y)) + geom_point(aes(colour = x1)) + xlab("Var.min") + ylab("Var.max") + ggtitle(title.main) +
scale_colour_gradient(low = "#00FFFF", high = "#000099")
#it becomes harder to seeplot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, -0.2),
plot.minmax.var = TRUE)
plot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, 0.2),
plot.minmax.var = TRUE)
re <- plot.ambcen.var(sd1 = c(0.2, 2), sd2 = c(1,1.8), rho = c(-0.5, 0.5),
weight.obj = 0.5, plot.minmax.var = FALSE,
return.ind = TRUE, plot.ind = TRUE)
#check the break point
weight.seq <- seq(0.1, 0.9, .1)
x.opt.seq <- numeric(length(weight.seq))
for (i in seq_along(weight.seq)){
w <- weight.seq[i]
re <- plot.ambcen.var(sd1 = c(0.2, 2), sd2 = c(1,1.8),
rho = c(-0.5, 0.5),
weight.obj = w, plot.minmax.var = FALSE,
return.ind = TRUE, plot.ind = FALSE)
x.opt.seq[i] <- re$x1.opt
}plot(weight.seq, x.opt.seq, type = "l") #it will go through a dramatic change
re <- plot.ambcen.var(sd1 = c(0.2, 2), sd2 = c(1,1.8), rho = c(-0.5, -0.4),
weight.obj = 0.5, plot.minmax.var = FALSE,
return.ind = TRUE, plot.ind = TRUE)
any(is.na(re$var.cen.seq))## [1] FALSEany(is.na(re$var.amb.seq))## [1] FALSEx.opt.seq2 <- numeric(length(weight.seq))
for (i in seq_along(weight.seq)){
w <- weight.seq[i]
re <- plot.ambcen.var(sd1 = c(0.5, 3), sd2 = c(1,3), rho = c(-0.5, -0.4),
weight.obj = w, plot.minmax.var = FALSE,
return.ind = TRUE, plot.ind = FALSE)
x.opt.seq2[i] <- re$x1.opt
}plot(weight.seq, x.opt.seq2, type = "l") #it is not a dramatic change 
#the optimal portfolio is robust to the choice of the weight in the objective function (the sharper or narrower the nose is, the less sensitive it will become.)6.4 Dynamic programming
#write a function by ourselves
h1 <- function(sig2, x1 = 0.5, rhol = -0.5,
sd1 = c(0.5, 2.0)){
re <- optimize(function(sig1) H(sig1=sig1, sig2=sig2,
rho = rhol, x1 = x1),
interval = sd1)
re$objective
}
h1(2)## [1] 0.75plot(Vectorize(h1), from = 0.1, to = 5)
#create a shiny app if needed